clear
set mem 1g
set more off
prog drop _all





** call program to generate moments 

global dir_data_out     "out_figs/" 


capture program drop getGroupMean
program getGroupMean

	args yvar xvar
	
	preserve
	
		keep if `yvar' != . & `xvar' != . 
				
		local ly : variable label `yvar'
		local lx : variable label `xvar'
		
		gcollapse (mean) `yvar', by(`xvar') 
					
		** output 
		outsheet `xvar' `yvar' using $dir_data_out/mean_`yvar'_BY_`xvar'.txt, nonames noquote nolabel replace
			
	restore
	
end



use out_files/sim_model_part1, clear
gen hgc = educ 

gen dq0 = dq if addiction == 0
label var dq0 "Initiate Smoking" 

gen logearnings_de0 = log(earnings) if enroll == 0 & earnings > 0 
label var logearnings_de0 "Log Earnings"


foreach yvar of varlist addiction dq0 hgc logearnings_de0 {
getGroupMean `yvar' age
}





******************************************************************
******************************************************************

** call program to generate figure 3


global dir_basedata    "../nlsy97/out_datasmm_new"
global dir_work_out    "out_figs"
global dir_fig_out     "../outputs"

global stats            "mean"
set scheme s1color



capture program drop plotfit 
program plotfit

	args yvar xvar


********************************************************************
***** Import first set of moments based on NLSY97 data

local stats = "$stats"

insheet using $dir_basedata/`stats'_`yvar'_BY_`xvar'.txt, clear

rename v1 `xvar'
*label var `xvar' "`lx'" 
*label values `xvar' _lxvalue

gen model = 1
rename v2 v_1

order `xvar' v_1 

save $dir_work_out/temp_data, replace

********************************************************************
***** Import second set of moments based on model simulated data

insheet using $dir_work_out/`stats'_`yvar'_BY_`xvar'.txt, clear

rename v1 `xvar'
*label var `xvar' "`lx'" 
*label values `xvar' _lxvalue

gen model = 2
rename v2 v_2

order `xvar' v_2


********************************************************************
***** Merge both datasets

merge 1:1 `xvar' using $dir_work_out/temp_data

*assert _merge == 1 | _merge == 3
keep if _merge == 3
drop _merge

save $dir_work_out/temp_data, replace




local label1 "NLSY Data"
local label2 "Benchmark Model"
local fig_se = 1


********************************************************************
***** Import first dataset

insheet using $dir_basedata/se`stats'_`yvar'_BY_`xvar'.txt, clear

rename v1 `xvar'
label var `xvar' "`lx'" 
label values `xvar' _lxvalue

gen model = 1
rename v2 semean

order `xvar' semean 


merge 1:1 `xvar' using $dir_work_out/temp_data 

*assert _merge == 1 | _merge == 3
keep if _merge == 3
drop _merge

gen v_1_min = v_1 - 1.96*semean 
gen v_1_max = v_1 + 1.96*semean


label variable `xvar' $lx 



**********************************************
*** Generate graphs

	local graph0_options "ms(i) c(l) lp(dot) lc(blue)"
	local graph1_options "ms(i) mc(blue) mfc(blue) c(l) lp(dash) lc(blue)"
*	local graph1_options "ms(T) mc(blue) mfc(blue) c(l) lp(dash) lc(blue)"
	local graph4_options "ms(T) mc(blue) mfc(none) c(l) lp(solid) lc(blue)"
	
	local graph2_options "ms(D) mc(red) mfc(red) c(l) lp(solid) lc(red)"
	local graph3_options "ms(O) mc(green) mfc(green) c(l) lp(solid) lc(green)"
	local graph5_options "ms(D) mc(red) mfc(none) c(l) lp(dash) lc(red)"
	local graph6_options "ms(O) mc(green) mfc(none) c(l) lp(dash) lc(green)"
	
		sum `xvar'
		local xmin = `r(min)'
		local xmax = `r(max)'
		local xsup = `xmax'+1 
		local xinf = `xmin'-1 		
	
		local step = 1 + 2*(`xmax'-`xmin'>20) 


	twoway rarea v_1_min v_1_max `xvar' , color(gs10) || scatter v_1 `xvar', `graph1_options' || scatter v_2 `xvar', `graph4_options' legend(lab(1 "95% CI" ) lab(2 `label1') lab(3 `label2')  ring(1)) ///
	ytitle(" $ly ") legend(on region(lstyle(none))) scheme(s1color) ///
	xlabel(`xmin'(`step')`xmax', valuelabel)  xscale( range(`xinf', `xsup')) ///
	ylabel($yinf($ystep)$ysup, grid) yscale( range($yinf, $ysup) ) 
	
    graph export $dir_work_out/fit_`stats'_`yvar'_BY_`xvar'.pdf, as(pdf) replace 
    graph export $dir_work_out/fit_`stats'_`yvar'_BY_`xvar'.eps, as(eps) replace 	
	

		
end 





** years smoked 
global yinf = 0
global ysup = 7
global ystep = 2
global ly "Yrs Smoked"
global lx "Age"
plotfit addiction age
graph export $dir_fig_out/Figure_03a.pdf, as (pdf) replace 
graph export $dir_fig_out/Figure_03a.eps, as (eps) replace 


** Figure 3.b
** smoke initiation rate 
global yinf = 0
global ysup = 0.15
global ystep = 0.05
global ly "Initiate Smoking"
global lx "Age"
plotfit dq0 age
graph export $dir_fig_out/Figure_03b.pdf, as (pdf) replace 
graph export $dir_fig_out/Figure_03b.eps, as (eps) replace 	


** Figure 3.c
** schooling over age
global yinf = 9 
global ysup = 16
global ystep = 1
global ly "Highest Grade Completed" 
global lx "Age"
plotfit hgc age
graph export $dir_fig_out/Figure_03c.pdf, as (pdf) replace 
graph export $dir_fig_out/Figure_03c.eps, as (eps) replace 	


** Figure 3.d
** log earnings when not in school 
global yinf = 8
global ysup = 12
global ystep = 1
global ly "Log Earnings"
global lx "Age"
plotfit logearnings_de0 age 
graph export $dir_fig_out/Figure_03d.pdf, as (pdf) replace 
graph export $dir_fig_out/Figure_03d.eps, as (eps) replace 	

!rm $dir_work_out/temp_data.dta 

capture log close 


